# Replication File for
# "Coalition Inclusion Probabilities: A Party-Strategic
# Measure for Predicting Policy and Politics
# Authors: Mark A. Kayser, Matthias Orlowski
# and Jochen Rehmert
# Code to replicate Table 5 and G1 and G2 (both Appendix)


# load required packages
# some of these may need to be installed first
library(stargazer);library(openxlsx)
library(multiwayvcov);library(Metrics)
library(car);library(lmtest);library(sandwich)

# set your directory to the folder with
# "data_deficit.RDS" data file 
setwd('...')

# load data
dat <- readRDS("data_deficit.RDS")

# prepare control variables
dat$minority <- 0
dat$minority[dat$cab_type == "minority"] <- 1

dat$periods <- NA
dat$periods[dat$year >= 1990 & dat$year <= 1994] <- 1
dat$periods[dat$year >= 1995 & dat$year <= 1999] <- 2
dat$periods[dat$year >= 2000 & dat$year <= 2004] <- 3
dat$periods[dat$year >= 2005 & dat$year <= 2009] <- 4
dat$periods[dat$year >= 2010 & dat$year <= 2014] <- 5
dat$periods[dat$year >= 2015 & dat$year <= 2019] <- 6


# Polls PM
summary(m.1 <- lm(deficit_diff ~ as.factor(country) +gdp_l + minority + yrs_next_election +num_pty + 
                    cab_range  + poll_pm_mean + periods,
                  data = dat[!is.na(dat$cip_gross_fin) & !is.na(dat$seatsh_fin),]))
m.1.se <- coeftest(m.1, vcov = vcovHC(m.1, "HC1"))[,2]

# Polls FIN
summary(m.2 <- lm(deficit_diff ~ as.factor(country) +gdp_l + minority + yrs_next_election +num_pty + 
                    cab_range  + poll_fin_mean +  periods  ,
                  data = dat[!is.na(dat$cip_gross_fin) & !is.na(dat$seatsh_fin),]))
m.2.se <- coeftest(m.2, vcov = vcovHC(m.2, "HC1"))[,2]

# Seats PM
summary(m.3 <- lm(deficit_diff ~ as.factor(country) +gdp_l + minority + yrs_next_election +num_pty + 
                    cab_range  + seatsh_pm + periods ,
                  data = dat[!is.na(dat$cip_gross_fin) & !is.na(dat$seatsh_fin),]))
m.3.se <- coeftest(m.3, vcov = vcovHC(m.3, "HC1"))[,2]

# Seats FIN
summary(m.4 <- lm(deficit_diff ~ as.factor(country) +gdp_l + minority + yrs_next_election +num_pty + 
                    cab_range +seatsh_fin + periods,
                  data = dat[!is.na(dat$cip_gross_fin) & !is.na(dat$seatsh_fin),]))
 
m.4.se <- coeftest(m.4, vcov = vcovHC(m.4, "HC1"))[,2]

 
# CIP PM
summary(m.5 <- lm(deficit_diff ~ as.factor(country) +gdp_l + minority +yrs_next_election +num_pty + 
                    cab_range  + cip_gross_pm + periods,
                  data = dat[!is.na(dat$cip_gross_fin) & !is.na(dat$seatsh_fin),]))
m.5.se <- coeftest(m.5, vcov = vcovHC(m.5, "HC1"))[,2]

# CIP FIN
summary(m.6 <- lm(deficit_diff ~ as.factor(country) +gdp_l + minority + yrs_next_election +num_pty + 
                    cab_range  + cip_gross_fin + periods,
                  data = dat[!is.na(dat$cip_gross_fin) & !is.na(dat$seatsh_fin),]))
m.6.se <- coeftest(m.6, vcov = vcovHC(m.6, "HC1"))[,2]

m.6$vcov <- cluster.vcov(m.6, ~ country)
m.6.se.rob <- coeftest(m.6, m.6$vcov)[,2]

# ----------------------------- #
# ----- Data for Table G1 ----- #
# ----------------------------- #
# vif
library(car)
mean(vif(m.2)[,3])
mean(vif(m.4)[,3])
mean(vif(m.6)[,3])

# RMSE
library(Metrics)
rmse(m.2$model$deficit_diff, m.2$fitted.values)
rmse(m.4$model$deficit_diff, m.4$fitted.values)
rmse(m.6$model$deficit_diff, m.6$fitted.values)

# heteroskedasticity
library(lmtest)
bptest(m.2)
bptest(m.4)
bptest(m.6)

# influential observations
mean(cooks.distance(m.2),na.rm = T)
mean(cooks.distance(m.4),na.rm = T)
mean(cooks.distance(m.6),na.rm = T)


# ------------------- #
# ----- Table 5 ----- #
# ------------------- #

stargazer(list(m.1, m.2, m.3, m.4, m.5, m.6),
          se = list(m.1.se, m.2.se, m.3.se, m.4.se, m.5.se,
                    m.6.se), omit = c("country", "periods"),
          covariate.labels = c("GDP Growth_{t-1}",
                               "Minority Cabinet",
                               "Time to Next Regular Election",
                               "Number of Cabinet Parties",
                               "Ideological Range in Cabinet",
                               "Polls PM (yearly average)",
                               "Polls FIN (yearly average)",
                               "Seatshare PM",
                               "Seatshare FIN",
                               "CIP PM (yearly average)",
                               "CIP FIN (yearly average)",
                               "Intercept"),
          out = "table5.tex")

###############################################################################
#           robustness checks
###############################################################################

# CIP PM
summary(m.7 <- lm(deficit_diff ~ as.factor(country) +gdp_l + minority +yrs_next_election +num_pty + 
                    cab_range  + cip_gross_pm + periods + keynes_cab_range,
                  data = dat[!is.na(dat$cip_gross_fin) & !is.na(dat$seatsh_fin),]))
m.7.se <- coeftest(m.7, vcov = vcovHC(m.7, "HC1"))[,2]


# CIP FIN
summary(m.8 <- lm(deficit_diff ~ as.factor(country) +gdp_l + minority + yrs_next_election +num_pty + 
                    cab_range  + cip_gross_fin + periods + keynes_cab_range,
                  data = dat[!is.na(dat$cip_gross_fin) & !is.na(dat$seatsh_fin),]))
m.8.se <- coeftest(m.8, vcov = vcovHC(m.8, "HC1"))[,2]


# CIP PM
summary(m.9 <- lm(deficit_diff ~ as.factor(country) +gdp_l + minority +yrs_next_election +num_pty + 
                    cab_range  + cip_gross_pm + periods + keynes_pm,
                  data = dat[!is.na(dat$cip_gross_fin) & !is.na(dat$seatsh_fin),]))
m.9.se <- coeftest(m.9, vcov = vcovHC(m.9, "HC1"))[,2]

# CIP FIN
summary(m.10 <- lm(deficit_diff ~ as.factor(country) +gdp_l + minority + yrs_next_election +num_pty + 
                     cab_range  + cip_gross_fin + periods + keynes_fin,
                   data = dat[!is.na(dat$cip_gross_fin) & !is.na(dat$seatsh_fin),]))
m.10.se <- coeftest(m.10, vcov = vcovHC(m.10, "HC1"))[,2]

# CIP PM
summary(m.11 <- lm(deficit_diff ~ as.factor(country) +gdp_l + minority +yrs_next_election +num_pty + 
                     cab_range  + cip_gross_pm + periods + rile_pm,
                   data = dat[!is.na(dat$cip_gross_fin) & !is.na(dat$seatsh_fin),]))
m.11.se <- coeftest(m.11, vcov = vcovHC(m.11, "HC1"))[,2]

# CIP FIN
summary(m.12 <- lm(deficit_diff ~ as.factor(country) +gdp_l + minority + yrs_next_election +num_pty + 
                     cab_range  + cip_gross_fin + periods + rile_fin,
                   data = dat[!is.na(dat$cip_gross_fin) & !is.na(dat$seatsh_fin),]))
m.12.se <- coeftest(m.12, vcov = vcovHC(m.12, "HC1"))[,2]



# -------------------- #
# ----- Table G2 ----- #
# -------------------- #

stargazer(list(m.7, m.8, m.9, m.10, m.11, m.12, m.6),
         se = list(m.7.se, m.8.se, m.9.se, m.10.se, m.11.se, m.12.se, m.6.se.rob), 
          omit = c("country", "periods"),
         add.lines = list(c("Number of Countries", length(unique(m.7[["model"]]$`as.factor(country)`)),
                            length(unique(m.8[["model"]]$`as.factor(country)`)),
                            length(unique(m.9[["model"]]$`as.factor(country)`)),
                            length(unique(m.10[["model"]]$`as.factor(country)`)),
                            length(unique(m.11[["model"]]$`as.factor(country)`)),
                            length(unique(m.12[["model"]]$`as.factor(country)`)),
                            length(unique(m.6[["model"]]$`as.factor(country)`)))),
          covariate.labels = c("GDP Growth_{t-1}",
                               "Minority Cabinet",
                               "Time to Next Regular Election",
                               "Number of Cabinet Parties",
                               "Ideological Range in Cabinet", 
                               "CIP PM (yearly average)",
                               "CIP FIN (yearly average)",
                               "Keynesian Policy: Range in Cabinet",
                               "Keynesian Policy: PM Position",
                               "Keynesian Policy: FIN Position",
                               "RiLe PM","RiLe FIN",
                               "Intercept" ),
         out = "tableG2.tex")
